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Abstract 

Several new discrete surface integral (DSI) methods for solving Maxwell’s equations in 
the time-domain are presented. These methods, which allow the use of general nonorthogo- 
nal mixed-polyhedral unstructured grids, are direct generalizations of the canonical staggered- 
grid finite difference method. These methods are conservative in that they locally preserve 
“divergence” or charge. Employing mixed polyhedral cells, (hexahedral, tetrahedral, etc.) 
these methods allow more accurate modeling of non-rectangular structures and objects be- 
cause the traditional “stair-stepped” boundary approximations associated with the orthog- 
onal grid based finite difference methods can be avoided. Numerical results demonstrating 
the accuracy of these new methods are presented. 
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Introduction 

The solution of physical problems whose behavior is governed by Maxwell's equations 
has been of considerable interest for many years. The propagation of electromagnetic 
signals such as microwaves for communication or radar pulses for the detection of aircraft 
are two examples of such problems. The conventional approach for numerically solving 
Maxwell’s equations in the time domain has been the use of finite difference methods 
(FDTD) in conjunction with orthogonal grids [1-5]. It is well known that the use of such 
methods can produce very accurate results (particularly when the domain is rectangular 
or almost rectangular). However, many problems involve domains or objects which are 
sufficiently irregular in their shapes, that they cannot be easily approximated with grids 
consisting of orthogonal cells [6]. 

Over the past years a number of approaches have been developed to treat problems 
with irregular shapes or boundaries. The simplest of these approaches has been to use 
the standard FDTD algorithm and to replace the curved and irregular boundaries with 
“stair-stepped” approximations. This approach has worked well in some cases but can give 
rather poor approximations in others unless very fine discretizatons are used [6]. Another 
approach has been to map the original irregular domain into a rectangular one and then 
apply the FDTD method in the transformed coordinates [7,8]. This approach suffers 
from the facts that more complicated equations must be solved and appropriate mapping 
functions are not always easily obtained. These approaches have been motivated by a 
desire to retain the use of the basic orthogonal grid FDTD method. The FDTD method 
has a number of attractive numerical properties. It is charge or divergence preserving 
both locally and globally. It is second order accurate and non-dissipative. It is also very 
computationally inexpensive. 

More recently, other numerical algorithms which allow the use of unstructured and 
irregular grids have been studied. Finite element methods, which have been very suc- 
cessfully used for solving for elliptic and parabolic partial differential equations, have been 
used with Maxwell’s equations [9,10]. These methods allow the use of boundary conforming 
nonorthogonal and unstructured grids which give one the ability to approximate irregu- 
lar boundaries with much greater geometrical accuracy. These methods usually provide 
for global charge or divergence conservation but not for local element charge conserva- 
tion. When applied on orthogonal structured grids, these methods do not reduce to be 
the familiar FDTD method. In fact, some of the finite element methods, when used with 
orthogonal grids, reduce to methods which have known difficulties such as grid decoupled 
solutions [9]. 

In a desire to use the FDTD algorithm essentially everywhere in a problem and still be 
able to match irregular boundaries, Taflove[ll] has developed an approximation technique 
that allows the use of special nonorthogonal cells immediately adjacent to the irregular 
boundary. Local conservation of charge for these special boundary cells is not preserved and 
second order accuracy is not assured in general. The primary advantage of this approach 
is that of very low computational cost. Long term stability of the method is also an open 
question. 

Drawing upon experience in computational fluid dynamics, finite volume methods have 
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recently been introduced to electromagnetics [12]. These methods use flux matching tech- 
niques across cell faces to form approximations to Maxwell’s equations. In order to achieve 
stability, these methods require the use of time integration methods which are dissipative. 
Lax-Wendroff or upwind time differencing methods have typically been used. These meth- 
ods preserve charge or divergence only to the truncation error level. Even though accurate 
short time solutions have been demonstrated, there remain accuracy questions. The longer 
term effects of the dissipative time integration methods and the effects of the lack of exact 
charge conservation are questions that remain to be answered. These methods also do not 
reduce to be the FDTD method when orthogonal grids are used. 

Not related to the above mentioned CFD based finite volume methods, but sharing a 
similar name, modified finite volume (MFV) methods have been recently presented [13j. 
These methods again allow the use of unstructured nonorthogonal boundary conforming 
grids and reduce to be the traditional FDTD method when used with structured orthogonal 
grids. These methods are charge conservative and non-dissipative and can produce very 
accurate solutions. Our first impression of these methods was that they had all of the 
desirable numerical approximation properties (charge conservation, non-dissipative, direct 
generalizations of FDTD, etc.). However, for some problems using highly distorted grids, 
weak instabilities and subsequent non-physical solution growth have been observed. 

With a continuing goal of developing algorithms without significant drawbacks for 
solving Maxwell’s equations using unstructured nonorthogonal grids, we seek algorithms 
which have the following properties: 

1. Allow the use of unstructured nonorthogonal grids. 

2. Allow a variety of cell or element types. 

3. Reduce to the FDTD method when orthogonal grids are used. 

4. Preserve charge or divergence locally (and globally). 

5. Conditionally stable. 

6. Non-dissipative. 

7. Accurate for nonorthogonal grids. 

In this paper we will present three new algorithms which meet the above criteria. 
These methods are derived using a discrete surface integration (DSI) technique. One might 
refer to these methods as finite volume methods, as they axe derived from an integration 
technique. However, since they involve only surface integrations, we prefer the terminology 
of discrete surface integration techniques. As formulated, the DSI techniques can be used 
with essentially arbitrary unstructured grids composed of convex polyhedral cells. ln ° ux 
implementation of the DSI algorithms, we allow for the use of unstructured grids whic 
are composed of combinations of nonorthogonal hexahedrons, tetrahedrons, triangular 
prisms and pyramids. These new methods actually reduce to be the conventional FD ID 
method when applied on an structured orthogonal hexahedral grid. The DSI techniques 
are formulated so that local (and hence global) conservation of charge or divergence of 
the fields is inherent in the algorithm and easily demonstrated. They use a leapfrog time 
integration technique which is conditionally stable and non-dissipative. This paper extends 
and improves our previous unstructured grid MFV algorithms [9,13]. The new algorithms 
have comparable accuracy and overcome the stability or solution growth difficulties which 
occurred when the MFV algorithms were applied on highly distorted grids. 
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It goes without saying that algorithms based on unstructured nonorthogonal grids 
will be significantly more complex and costly than the orthogonal grid FDTD method. 
However, since the DSI methods reduce to be the computationally efficient FDTD method 
when the grid is orthogonal, we advocate the use of grids which are orthogonal almost 
everywhere so that advantage may be taken of the simplifications that occur in this case. 
Even when nonorthogonahty exists, the use of efficient data structures allow the new DSI 
algorithms to be quite computationally efficient. 

In the next section, we will describe in detail the new DSI methods. We will then 
discuss issues concerning the implementation of these methods as they relate to efficiency. 
We will then describe the numerical experiements that have been performed to assess 
stability. Finally, numerical results will be presented to demonstrate the accuracy of the 
DSI methods. 


Discrete Surface Integration Methods 

We begin by assuming that we wish to solve Maxwell’s curl equations on an irregular 
three-dimensional domain R which has a boundary surface denoted by S. We will also 
assume that the domain R has been discretized into convex polyhedrons. Fig. 1 shows a 
twisted waveguide discretized using hexahedral cells. This is an example of a problem type 
that we wish to consider that could not be easily solved using the conventional orthogonal 
grid FDTD method. Maxwell’s curl equations are given by: 


9D 

dt 

dB 

dt 


V x H 

(1) 

-V x E 

(2) 


where for linear isotropic materials the vectors D, E, B and H are related by the consti- 
tutive relationships 


D = cE 

B = /iH. 

The linear isotropic material properties are: e, the permittivity, and /z, the permeability. 
Like the conventional FDTD method, the DSI methods can be generalized to treat more 
complex materials. However, for the purposes of this paper, we will assume that the 
linear isotropic material properties are piecewise constant over the domain R. We will 
also assume that 5 is a perfectly conducting surface, i.e., Eta , n = 0. This is sufficient to 
guarantee that the problem is well-posed. 

We will derive our new DSI algorithms for unstructured grids formed from convex 
polyhedral cells. We restrict the choice of cell types to convex polyhedrons whose edges 
are straight lines. The faces of the polyhedrons are not necessarily planar and we make 
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the assumption that any face in the assembled grid is shared by at most two cells. These 
are very weak restrictions and allow a great deal of flexibility. 

The DSI method requires the use of a dual grid. The dual grid and its structure are 
completely derivable from a knowledge of the primary grid. For each cell of the primary 
grid, we define its barycenter to be a node of the dual grid. The barycenter of a cell is 
located at the average of the coordinates of the nodes which define the cell. We construct 
edges of the dual grid by connecting barycenters of adjacent cells with straight lines passing 
through each of the interior cell faces of the primary grid. The barycenters of two cells will 
be connected (i.e., form a dual edge) if and only if the two cells share a common face. For 
primary cell faces which lie on the problem boundary, S , we form a corredponding dual edge 
by joining the cell barycenter to the face barycenter. There is a one-to-one correspondence 
between the nodes, edges, faces, and cells of the primary grid to the cells, faces, edges and 
nodes of the dual grid, respectively. The dual face associated with a primary edge has as 
its perimeter the dual edges associated with all of the primary faces which share the given 
primary edge. The dual cell associated with a primary node has as its surface the dual 
faces which correspond to all of the primary edges which share the primary node. Though 
not necessary for the definition of the new algorithms, we recommend that the variations 
in grid sizes and angles be sufficiently smooth so that the primary and dual edges actually 
intersect their corresponding dual and primary faces, respectively. Degradation of solution 
accuracy has been observed when this condition is not met. Fig. 2 shows an eight cell 
hexahedral primary grid and its one interior dual cell. 

Our DSI solution variables will be associated with the edges and faces of the primary 
grid and also with the edges and faces of the dual grid. The quantity associated with a 
pr im ary cell edge is the projection of the electric field vector onto that edge, i.e., E-s, where 
s is the primary cell edge vector. The magnetic field projection H • s* is associated with 
a dual cell edge where s* is the dual cell edge vector. In addition, with each primary grid 
face we will associate a full magnetic field vector B, and with each dual grid face we will 
associate a f ull electric displacement vector D. We will denote with an asterisk, *, geometric 
quantities associated with the dual grid. Fig. 3 depicts these associations. We remark that 
these associations of field quantities with the primary and dual grid locations are entirely 
reciprocal and that the respective locations of the magnetic and electrtic field quantities 
could be interchanged. The particular choice of which field quantities to associate with 
each grid is best determined by deciding which field quantities one desires to have on the 
exterior boundary surfaces where the boundary conditions will be imposed. Since we are 
assuming that our domain R is surrounded by a perfect electric conductor, we will associate 
the electric and magnetic field quantities as described above. For open region problems, 
the choice for the location of the field quantities will depend on the particular radiation 
boundary condition algorithm used. 

We will now describe the equations and algorithmic process used to advance in time the 
magnetic field vectors which are associated with each primary grid face. We assume that 
the time variable has been discretized by the choice of a timestep size, At. Superscripts 
on field quantities will denote their time state with D* = D(<*) = D(fcAf). As we 
will be using a leapfrog style time integration method, the magnetic field vectors will be 
associated with half-integer times, t k+ >, and the electric field vectors will be associated 
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with integer times, t k . For a particular primary cell face, we define the area-normal vector 
to be N = / n dS, where n is a unit surface normal defined by the right-hand rule in 
relation to a specified circulation around the perimeter of the cell face. It is easily shown 
that N is uniquely determined by the perimeter of the surface and is independent of the 
actual interior surface shape. This fact allows for simple computations of N using piecewise 
planar approximating surfaces. Using (2) for the primary grid face, F , we define the time 
derivative of the normal component of the magnetic field to be: 



( 3 ) 


Equation (3) allows us to obtain the time derivative of the normal component of the 
magnetic field on a primary face, F, from the electric field projections onto the perimeter 
edges of that face. The last integral in (3) is easily computed numerically by summing 
these edge projections. This can be done for each primary cell face. 

The next step in the algorithm is to use these time derivatives of the normal components 
of the magnetic field to compute a full vector value of dBjp / dt for the primary cell face, 
F. We will assume that the face F is defined by P primary edges and nodes, with the i th 
node being located at the intersection of the consecutive edges z and m = (z mod P) + 1. 
Also, we assume that the face F is shared by N c primary grid cells, where by construction: 
N c = 1 for boundary faces and N c = 2 for interior faces. We will denote by Fij , the face 
of cell j (other than F) which shares edge z. Fig. 4 depicts these associations for a dual 
edge associated with a primary face defined by five primary edges. At each of the P nodes 
of face F, we will unfold N c vector values <£B k j/dt by solving the 3x3 system of equations: 


dt 


N f = 


d B 


dt 

dt 


M*, = 


N F m , = 


- f E k • dl 

- / E k • dl 

-l E k • dl , 
J F* n..i 


( 4 ) 


where z = 1, . . . , P, j = 1, . . . , N e , and m = (i mod P) + 1. For this primary face, F , which 
is shared by N c primary cells, there are PN C different values of dB^/dt, which will now 

be averaged or interpolated to form a single dB k F /dt vector for the face. There are many 
reasonable ways to average these vectors and we will consider the three following averaging 
or interpolating methods: 


dBj? 

dt 


= ( 


N C P 


dB fc , 

>££ — 

j = 1 1=1 


dt 


(5a) 
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dB k F 

dt 


dt 


dB k F 

dt 



N c P 

S £ to I 

; = li=l 


dt 


N c P 

EE to I 

j=li=l 


? 


( 56 ) 


(5c) 


where the weight 

wij = N p • (N F itj x N Fi + i,i) » 

represents the volume of the j** local coordinate system at node i of face F. We note 
also that the weight Wij is the determinant of the system of equations (4). A different dis- 
cretization method will be obtained for each different averaging method. We characterize 
(5a) as a simple vector sum average, (5b) as a partially volume weighted average, and (5c) 
as a fully volume weighted average. As will be subsequently discussed, numerical evidence 
indicates a preference for the fully volume weighted averaging scheme (5c). 

The full B vector for each primary face, F, may now be advanced in time using the 
time-centered leapfrog algorithm: 


B* + 5 = B*"5 + A , ( 6 ) 

dt 

where At is the specified timestep size. Finally, a time advanced value of the projection 
of the magnetic field onto the dual edge, s*, which penetrates the primary cell face, F, is 
easily obtained using: 


(H-s*) fc+ 5 


B fc+ i s* 

“ “ ? 

P- 


( 7 ) 


where n is an appropriate permeability value. If the permeability is discontinuous then a 
value can be determined by an appropriate average as is done for FDTD algorithms. 

It is important to note that the first equation in (4) (namely, the equation coming from 
the given face F) is common to all of the sets of 3 x 3 equations being used in the time 
advance process for the face. This implies that the averaged value (5) of the time derivative 
of the magnetic field vector for F also satisfies this equation. It is from this aspect of the 
numerical algorithm that divergence or charge conservation can be demonstrated. If we 
integrate the divergence of (6) over a primary cell, C , we obtain 
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J(V-B k+ i)dV = j(V- 

Considering the last integral term in (8) we have 


(8) 


dV + At [ (V 

Jc 


dB 

dt 


k 

)d.V. 


Jv f» v - Jj ¥ " )ds - - - °' (9) 

where the sum runs over all the faces Fi of the closed cell C . The last sum is zero because 
each edge of C will be traversed twice, once in each direction. Thus we see that the local 
divergence of the time derivative of the magnetic field vector is zero and so if the initial 
fields have zero local divergence, then (8) and (9) prove that zero divergence of the time 
advanced fields will be conserved. 

To time advance the electric field vectors which are associated with each dual cell face, 
we proceed in a manner which is exactly “dual” to the magnetic field procedure described 
above. For any dual cell face, we define the dual area-normal vector to be N* — J n* dS% 
where n* is a unit dual surface normal defined by the right-hand rule in relation to a 
specified circulation around the perimeter of the dual cell face. Using (1) for a given 
dual grid face F *, we define the time derivative of the normal component of the electric 
displacement vector to be: 


d D t+ 5 
dt 


■ Np. 



( 10 ) 


Again, the last integral in (10) is easily computed by summing the projections of the 
magnetic field on the edges defining the dual cell face. This can be done for each dual cell 
face. 

The next step in the algorithm is to use these time derivatives of the normal components 
of the electric field to compute a full vector value of dDp./dt for the dual cell face, F* . We 
will assume that the face F* is defined by P * dual edges and nodes, with the i th dual node 
being located at the intersection of dual edges i and m = (t mod P*) + 1. Also, we assume 
that the dual face F* is shared by N* dual cells. We will denote by F*j , the dual face of 
dual cell j (other than F*) which shares dual edge t. Fig. 4 with all quantities replaced 
by their appropriate duals would depict these associations. At each of the P* dual nodes 

of dual face P*, we will unfold N* vector values dD^'/dt by solving the 3 x 3 system of 
equations: 
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<®‘ +! 

■ Nr. = 

I H k+ *dl* 


dt 

£ 

Jf * 


<iD* +S 

dt 

■ Np. = 

F i,j 

<f H k+ 5 • d\* 

Jf* 

(ID 





t; 

Nr. = 

I H k+ 5 • dl* , 


dt 

*">i) 

Jf* 

m,j 



where i = 1, . . . ,P*, j = 1,. . . , 1V C * and m = (i mod P*) + 1. For this dual face, F* , which 

is shared by N* dual cells, there are P*N c * different values of dD-/ 5 /^» which wiU now be 

averaged or interpolated to form a single dD p.^ / dt vector for the dual face. The averaging 
methods (5a-5c) used for the magnetic fields are also used for the electric fields. 

The full D vector for the dual face, F*, may now be advanced in time using the 
time-centered leapfrog algorithm: 


jrk fc+i 

D*+l = D* + A/-r- , ( 12 ) 

at 

where At is the specified timestep size. Finally, a time advanced value of the projection of 
the electric field onto the primary edge, s, which penetrates the dual cell face, F*, is easily 
obtained using: 


(E • s ) fc+1 


D * +1 • s 

€ 


(13) 


where e is an appropriate permittivity value. The local conservation of divergence for 
the electric field can be easily shown in a manner similar to that described above for the 
magnetic fields. 

Equations (3)-(7) for the magnetic field quantities, equations (10)-(13) for the electric 
field quantities and a linear averaging method (5), constitute the new divergence conserving 
DSI approximation methods. 

When used with structured grids composed solely of hexahedral cells, the new DSI 
method differs only slightly from our previous MFV method. In computing the time 
advanced value of an edge-projected magnetic field value, exactly the same electric field 
values would be used to compute the required time derivative. However, slightly different 
weights or coefficients would be used in the DSI method. For unstructured grids, the 
differences between the DSI method and the MFV method become more pronounced. 
In general, the DSI method requires fewer dual edge magnetic field projections to time 
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advance a given electric field projection than does the MFV method. For the DSI method, 
only those dual edges defining the dual faces which are immediately adjacent to the dual 
face of the primary grid edge being updated are used in the update algorithm. In contrast, 
the MFV algorithm would use the dual edges of all of the dual faces (potentially many 
more) of the two dual cells which surround the primary grid edge field value that is being 
time advanced. Thus we see that the DSI algorithm is in general more compact than is 
the MFV algorithm. 

We note that significant simplifications occur when a primary edge is orthogonal to 
its dual face, or “dually”, when a dual edge is orthogonal to its primary face. When this 
occurs, the vectors s and NJ-. (or s* and N f) are aligned. The time advance of the edge 
projected field value may be performed directly using (3) (or (10)) and the averaging (5) 
may be completely bypassed. Stated another way, when the above orthgonality conditions 
exist, the time advance of E • s (or H • s*) may be accomplished using a single line integral 
around the perimeter of the face F * (or F). It is this fact that demonstrates that the 
DSI methods are completely equivalent to the canonical FDTD methods when orthogonal 
hexahedral based grids are used. Therefore, the DSI methods are direct generalizations of 
the orthogonal grid FDTD methods for unstructured nonorthogonal grids. These orthog- 
onality conditions may occur locally for only a relatively few edges or may occur globally 
for the entire grid. They occur globally when structured grids composed of orthogonal 
hexahedral cells are used. They also occur globally when grids are used which are three 
dimensional analogs of two-dimensional grids formed from Delaunay triangles and their 
dual Voronoi polygons. We will indicate in the next section how a properly structured 

code can take advantage of these simplifications when orthogonality occurs. 

• 

Exterior (perfect electric conductor) boundary conditions are very easily handled by 
simply setting E • s = E tan = 0 for those edges s of the grid which form the boundary. 
Therefore, the above time advance procedure is not required for the electric field projections 
for edges which are boundary edges as these projections are determined completely from 
the boundary conditions. For primary cell edges which have one endpoint on the problem 
boundary surface 5 and the other in the interior of R, there will be only one complete 
dual cell associated with this edge (rather than two). The second dual cell is incomplete 
as it is clipped by the boundary surface 5. For computing E s for such edges, appropriate 
values for the fine integrals around the incomplete dual faces can be determined from 
the boundary conditions using reflection principles. By choosing primary grids which are 
locally orthogonal near the exterior boundary, the time advance of the edge projection 
ran be performed as mentioned above without using any of the incomplete boundary dual 
faces. 

Local grid orthogonality at exterior boundaries also facilitates the solution of open 
region problems (those which are not enclosed by perfect conductors). We have not as 
yet developed three-dimensional radiation boundary conditions which are useable with 
general nonorthogonal grids. However, open region problems can be solved using the new 
DSI methods if the exterior boundary cells are made to be orthogonal. In this case, the 
radiation boundary condition methods [14] developed for the canonical FDTD methods 
are easily applied. Again, this is possible because of the equivalency of the DSI methods 
and the canonical FDTD method for orthogonal grids. 
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As implied in the DSI algorithm description above, interior material interfaces are also 
easily handled by using appropriately averaged values for the discontinuous permittivity, 
c, and the permeability, fi, (as is customary with the conventional orthogonal grid FDTD 

method). 


Implementation of the DSI Method 

In actually implementing the new DSI algorithm, there is considerable flexibility in 
choosing the variables to be used as the actual computational variables. If computer stor- 
age were not a significant consideration, one could keep the edge projected values, the 
face normal components, and the full solution vectors associated with the faces. However, 
storage is almost always an issue for large three-dimensional problems and in our imple- 
mentation, we have chosen to use the projections of the vector fields onto the primary and 
dual grid edges. Direct equations for the time advance of H • s and E • s are easily derived 
by combining and manipulating (3)-(7) and (10)-(13), respectively. For an implementation 
using edge projected field values as the variables, the total number of unknown approxi- 
mate field quantities is clearly equal to the number of cell edges in the primary grid plus 
the number of cell edges in the dual grid. For a nonorthogonal but structured hexahedral 
grid, the time advance of a typical interior E • s will involve the use of 20 surrounding 
magnetic field projected values, H • s. A similar number of electric field values would be 
used for the time advance of a typical magnetic field value. For orthogonal hexahedral 
grids, this number reduces to 4. These significant potential savings are the basis for our 
advocating the use of orthogonal grids wherever possible. These potential savings have 
also motivated us to design an implementation of the algorithm which can benefit from 
orthogonality when it occurs. 

In implementing this new DSI algorithm, we have structured the code so that there first 
occurs a pre-processing phase. In this phase, the dual grid is constructed from the primary 
grid information and a dependency graph is constructed for each primary grid and dual 
grid variable. A dependency graph is simply a list of pointers which point to all of the other 
grid variables that are required to time advance a specific grid variable. In conjunction 
with the dependency graph, coefficients are also generated so that any given variable may 
be advanced in time by computing a single very simple dot product of a coefficient vector 
with a vector of appropriate variables determined from the dependency graph. In forming 
the dependency graph, checks for orthogonality are performed and variables that are not 
needed for the time advance due to reasons of orthogonality are never used or stored in 
the dependency graph. This allows us to take advantage of the greater efficiencies (both in 
storage and computational arithmetic) inherent in the algorithm when grid orthogonality 
exists, either locally or globally. Following the pre-processing phase is the execution phase 
where the dependency graph, coefficients, and simple dot products are repeatedly used to 
time advance all of the problem variables. 

We have recently been solving some rather large problems using the DSI method on 
unstructured grids. One particular problem involved an unstructured grid consisting of 
262,620 cells, resulting in 1,645,195 degrees of freedom, and required only about 17,493,000 
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words of memory to execute on a Cray-2. No memory saving techniques such as word pack- 
ing were used. The computational time required for this problem was about 7 microseconds 
per timestep per degree of freedom. The amount of memory and computational time re- 
quired for a problem can vary significantly depending upon the structure and orthogonality 
present in the problem grid. 


Stability Experiments 

We have presented a description of three new DSI algorithms for solving Maxwell’s 
equations. Each algorithm is distinguished by a different averaging method used to form 
the time derivative of a vector field associated with a cell face. As these methods reduce to 
be the canonical FDTD method on orthogonal grids, their stability properties (when used 
with this type of grid) are the same as for the FDTD methods. As we have been unable to 
mathematically analyze their stability properties for nonorthogonal unstructured grids, we 
have performed computational experiments to study this issue. While we recognize that 
these experiments do not conclusively demonstrate stability, we feel that they are useful 
to give general indications of stability characteristics. 

For testing stability on highly distorted grids, we have found that a grid obtained from 
mapping a 10 x 10 x 10 cubical grid to a sphere provides a good test. The corner cells 
of the cubical grid become very highly distorted under the mapping. For two dimensional 
testing, a square grid mapped to a circle provides a good test. Fig. 5 shows a 20 x 20 
two dimensional square grid mapped to a circle and indicates the degree of distortion that 
occurs in the corner cells. These problems can be randomly excited and the solution can 
be monitored for many thousands of time steps to see if growth occurs in the solution. We 
judge the stability of an algorithm by measuring the solution growth over time in terms of 
total energy or maximum solution magnitude. 

We have experimented with the three different DSI averaging methods: (5a), a simple 
vector sum average of all of the time derivative vectors on both sides of the given face; 
(5b), a vector sum average of the two volume weighted averages of time derivative vectors 
on either side of the given face; and (5c), a full volume weighted average of all of the time 
derivative vectors on both sides of the given face. Numerical experiments involving the 
use of 50000 time steps have indicated that all of the three DSI methods are conditionally 
stable. However, they differ in that they seem to have different stability limits. For 
the highly distorted spherical grid mentioned above, we have observed differences in the 
maximum allowable time step sizes by factors as large as five. The simple vector sum 
average method, (5a), always required the smallest time step size and the full volume 
weighted averaging method, (5c), allowed the largest time step size. The mixed vector sum 
and volume weighted average method, (5b), always fell between the other two methods. 
For these stability reasons we advocate the use of the full volume weighting averaging 
method (5c). In contrast, our previous MFV method[13] is simply unstable for these 
highly distorted grids. In performing stability experiments with various algorithms, we 
have observed that it is essential to use a large simulation time (involving many time steps) 
so as to allow enough time for signals to transit the computational domain many times. 
For some distorted grids, our previous MFV method would exhibit no discernible solution 
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growth until many thousands of time steps had occurred. Fig. 6 shows an example of an 
MFV solved waveguide problem that showed no apparent solution growth until some 20000 
time steps had been taken. In Fig. 6, the electric field at every tenth time step for a single 
point is plotted versus the time step number. For this same problem, the DSI algorithm 
does not show any spurious solution growth even to 50000 time steps for. To determine an 
appropriate time step size in solving a particular unstructured grid problem, we compute 
the smallest primary or dual grid edge length and divide by the local propagation velocity. 
We then use a time step size which is half of this computed value. In general, this has 
been a reliable approach. 


Accuracy Experiments 

We feel it is important to demonstrate the accuracy of new algorithms by compar- 
ing their numerical solutions for problems which have known and computable analytic 
solutions. We have extensively compared the new DSI method with our previous MFV 
method[13]. For problems using hexahedral grids where the MFV method does not en- 
counter stability related solution growth, we have found that the two methods produce 
solutions which are almost indistinguishable. Therefore, in the results which follow, we 
will present results only for the full volume weighted averaging DSI method. 

We first consider a simple problem which is designed to provide some assessment of the 
wave propagation characteristics of the new algorithm. We consider the propagation of a 
TEio signal in a rectangular waveguide which is 1 meter in width, 0.3 meters in height, and 
10.0 meters in length. The signal is initiated in the waveguide by specifying the tangential 
field, E z , at the left boundary (y = 0.0). The function used to drive this signal is: 


E z (i,x,0,z) = «tn(1.5iri) sin(irx) . (14) 

We assume that all of the waveguide walls are perfect electric conductors, i.e. Etan = 0. 
Normalized values for the material parameters e and n are used: e = 1.0, and n = 1.0. 
The signal is allowed to propagate until i = 8. The analytic solution for this problem has 
been derived and computed so we may examine the errors precisely. 

We have used a variety of grids to numerically solve this problem. Clearly one simple 
grid for this problem is an orthogonal hexahedral grid. However, to demonstrate the 
ability to use nonorthogonal grids, we deliberately skew the 10 x 100 x 3 grid so that 
it is composed of nonorthogonal hexahedrons as shown in Fig. 7. The skewing for this 
grid is such that halfway down the length of the guide the cells have angles of about 
57 and 123 degrees. Fig. 8 compares the time histories of the numerical solutions for 
the skewed and orthogonal hexahedral grids, respectively, with the analytic solution for a 
point in the center of the waveguide and 4.0 meters down the guide. We notice that the 
significant error occurs as the signal first reaches the observation point. Both solutions 
exhibit typical dispersion errors; however, the skewed grid results are somewhat better. 
To more precisely compare the solutions, Fig. 9 shows the maximum errors over the entire 
waveguide for both the skewed and orthogonal hexahedral grid solutions as a function 
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of time. These maximum errors should be referenced to the solution magnitude which 
oscillates roughly between +1 and —1. Surprisingly, for all but the earliest times, the 
skewed grid produces a uniformly better solution. To demonstrate that other grid cell 
types can produce acceptable solutions, Table I compares the maximum errors, run times 
in seconds for an IBM RISC 6000 Model 550 computer, and problem sizes for several 
different grids. The grids differ in the number of cells and edges. The non-hexahedral 
grids were all formed from the orthogonal and skewed hexahedral grids by subdividing 
each hexahedral cell into 5 tetrahedrons, 6 pyramids, 2 tetrahedrons and 2 pyramids, or 2 
triangular prisms, The subdivision of a single hexahedral cell into 5 tetrahedrons and the 
subdivision into 2 tetrahedrons and 2 pyramids are shown in Figs. 10-11. 

From Table I one can see that the most accurate results were obtained using the tetra- 
hedral grid derived from the skewed hexahedral grid. However, these results required a 
great deal more computer time to obtain. The worst results resulted from the use of the 
skewed grid consisting of tetrahedral and pyramidal cells. The poorer results reflect the 
fact the grid is composed of highly distorted cells which have very large and very small 
internal angles. It seems apparent that this problem is most efficiently solved using grids 
consisting of hexahedral cells. 

To demonstrate the accuracy that can be gained by the better zoning of nonorthogonal 
surfaces than can be achieved using only orthogonal cells, we consider the problem of the 
scattering of a plane wave gaussian pulse from a perfectly conducting sphere of radius 0.1 
meters. The plane-wave pulse is defined by the Gaussian-like function 



if 0 < t < 2 x 10“ 9 , 
otherwise. 


This problem again has a computable analytic solution so we may compare errors precisely. 
Fig. 12 shows the DSI surface grid and a stair-stepped orthogonal grid approximation 
for the surface of a sphere of radius 0.1 meters. Both surface grids are derived from 
a 6 x 6 x 6 cubical grid. The DSI surface grid is obtained by mapping the surface of 
the discretized cube to a 0.1 meter radius sphere. The stair-stepped orthogonal grid 
is obtained by elimin ating cells from the cubical grid. Fig. 13 compares the orthogonal 
grid finite difference scattered field solution and the DSI scattered field solution with the 
analytic solution for a point which lies in the shadow region (behind the sphere relative 
to the direction of the incident pulse) two sphere radii away from the sphere surface. The 
maximum errors are noticeably different. Significantly more accurate results are obtained 
by avoiding the orthogonal ‘‘stair-stepped” grids. It is also noteworthy that the DSI results 
were obtained using a 2.2 meter diameter hemispherical domain which was discretized using 
only 9600 radially expanding cells. In contrast, the finite difference solution used half of a 
2.0 meter cubical domain with a 60 x 60 x 30 uniform cubical discretization (108000 cells). 

Another problem that has a known analytic solution is that of a pulsed dipole radiating 
in free space. We assume that a dipole is located at the origin and is oriented in the z- 
direction. We will drive it for a finite amount of time with the pulsed amplitude function 
f(i) = p(a(f)), where a(<) = 4 x 10 7 - 1, and where 
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if -1 < X < 1, 
otherwise. 



As the solution of this problem depends strongly on this function and its first two deriva- 
tives, we have chosen the function to have smooth second derivatives. The analytic solution 
may be found in almost any standard text or in [15]. To solve this problem, we use a grid 
that is mostly orthogonal but in the corner nearest the origin we use a small volume of 
nonorthogonal cells which transition from a spherical surface to the cubical grid. Fig. 14 
shows the grid used. We excite the problem by specifying the tangential components of the 
analytic solution on the spherical boundary surface and compare the numerical solution 
with the analytic solution at a point just outside of the nonorthogonal region. Fig. 15 
compares E z component of the DSI solution with the analytic solution for a point located 
slightly outside of the nonorthogonal cell region. We see that excellent results are obtained. 
The other solution components exhibit comparable accuracy. 

As a final example we consider the problem of computing a TE \ o waveguide solution for 
the rectangular waveguide (shown in Fig. 1) which goes through a 180 degree twist. This 
problem cannot be effectively solved using a stair-stepped orthogonal grid finite difference 
algorithm unless a tremendously large number of cells are used. The analytic solution for 
this problem is not known. However, from actual waveguide experiments, we know that 
a TE\q mode should propagate and be rotated 180 degrees with little distortion by such 
a waveguide. Fig. 1 shows the 10 x 100 x 3 cell DSI grid used and Fig. 16 compares the 
DSI solution time history for a point in the center of the waveguide immediately after the 
twist with a TE \ o solution for a similarly discretized waveguide with no twist. We notice 
that the waves transition the twist with little distortion. Conversion of the TE\q mode to 
other modes does not occur because of the thinness of the waveguide. This problem was 
excited using the driving function (14). 


Conclusions 

The numerical results presented here lead us to believe that the DSI method pre- 
sented in this paper is one with significant advantages. First, it is a direct generalization 
for unstructured nonorthogonal grids of the well-known and widely used orthogonal grid 
FDTD method for Maxwell’s curl equations, so its approximation properties on regular 
orthogonal grids are quite well understood. Second, the DSI method is conditionally stable 
even for highly distorted grids. This is a distinct improvement when compared with our 
previous MFV algorithm. Third, the DSI method is conservative in that it conserves the 
divergence of the fields both locally and globally. Fourth, it is quite accurate (comparing 
well with our previous MFV algorithm) even when used with irregularly structured grids. 
Fifth, it is very flexible and can be used with various combinations of different polyhedral 
discretizations . 

The most significant disadvantage of the DSI method is the larger amount of informa- 
tion (primary grid structure plus dual grid structure) that is required. This disadvantage 
can largely be overcome by using discrete grids which are primarily composed of orthogonal 
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hexahedral cells together with a relatively few number of irregularly shaped nonorthogonal 
cells which are used to match boundaries more accurately. Implementing the algorithm in 
a two-phase way where all of the geometry information is processed first and condensed 
into a dependency graph and related coefficients allows one to take advantage of grid 
orthogonality when it occurs and also mitigates the extra nonorthogonal grid costs. 
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Figure Captions 


Figure 1. Twisted waveguide discretized into hexahedral cells. 

Figure 2. Primary grid consisting of eight hexahedral cells and its one interior dual cell. 

Figure 3. Discrete electric and magnetic field variable locations relative to the primary and dual 
grid cells. 

Figure 4. Primary grid faces used to time advance a magnetic field vector. 

Figure 5. Highly distorted grid obtained by mapping a 20 x 20 square grid to a circle. 

Figure 6. Long-term solution growth plot for the weakly unstable MFV algorithm. 

Figure 7. Skewed 10 x 3 x 100 nonorthogonal hexahedral waveguide grid. 

Figure 8. Numerical TE\q waveguide solutions compared with the analytic solution for a point 
4.0 meters down the guide: a) analytic solution, b) orthogonal grid DSI solution, c) 
skewed grid DSI solution. 

Figure 9. Comparison of maximum errors as a function of time for the DSI orthogonal and skewed 
waveguide hexahedral grids. 

Figure 10. Hexahedral cell subdivided into five tetrahedrons. 

Figure 11. Hexahedral cell subdivided into two tetrahedrons and two pyramids. 

Figure 12. Surface discretizations of a 0.1 meter radius perfectly conducting sphere: a) Boundary 
conforming DSI surface grid (left), b) stair-stepped orthogonal FDTD grid (right). 

Figure 13. Comparison of the DSI numerical solution and the stair-stepped orthogonal grid FDTD 
solution with the analytic solution for a point in the shadow region two sphere radii 
from the scatterer. 

Figure 14. Mostly orthogonal DSI grid for the dipole problem with a nonorthogonal region tran- 
sitioning from a spherical surface. 

Figure 15. Comparison of the E z component of the DSI and analytic dipole solutions for a point 
just outside the nonorthogonal grid region. 

Figure 16. Comparison of the electric field time histories for the twisted waveguide DSI numerical 
TEio solution with a non-twisted waveguide DSI numerical solution at a point in the 
center of the guide immediately beyond the twisted portion. 
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Table I. Comparison of TE \ o waveguide solution maximum errors, run times, and problem 
sizes for various DSI grids. 
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Figure 2. Primary grid consisting of eight hexahedral cells and its one interior dual cell. 
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Figure 4. Primary grid faces used to time advance a magnetic field vector. 
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Figure 5. Highly distorted grid obtained by mapping 
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Figure 7. Skewed 10 x 3 x 100 nonorthogonal hexahedral waveguide grid. 
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Figure 10. Hexahedral cell subdivided into five tetrahedrons. 



Figure 11. Hexahedral cell subdivided into two tetrahedrons and two pyramids. 




Figure 12. Surface discretizations of a 0.1 meter radius perfectly conducting sphere: a) Boundary 
conforming DSI surface grid (left), b) stair-stepped orthogonal FDTD grid (right). 
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Figure 14 . Mostly orthogonal DSI grid for the dipole problem with a nonorthogonal region tran- 
sitioning from a spherical surface. 
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Table I. Comparison of TE\o waveguide solution maximum errors, run times, and problem 
sizes for various DSI grids. 
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